Laws of crack motion and phase-field models 

of fracture 



Vincent Hakim ^ 

^Lahoratoire de Physique Statistique, CNRS-UMR8550 associe aux universites 
Paris VI et VII, Ecole Normale Superieure, 24 rue Lhomond, 75231 Paris, France 

Alain Karma ^ 

^ Physics Department and Center for Interdisciplinary Research on Complex 
Systems, Northeastern University, Boston, Massachusetts 02115 



Abstract 

Recently proposed phase- field models offer self-consistent descriptions of brittle frac- 
ture. Here, we analyze these theories in the quasistatic regime of crack propagation. 
We show how to derive the laws of crack motion either by using solvability condi- 
tions in a perturbative treatment for slight departure from the Griffith threshold, 
or by generalizing the Eshelby tensor to phase-field models. The analysis provides 
a simple physical interpretation of the second component of the classic Eshelby 
integral in the limit of vanishing crack propagation velocity: it gives the elastic 
torque on the crack tip that is needed to balance the Herring torque arising from 
the anisotropic interface energy. This force balance condition reduces in this limit 
to the principle of local symmetry in isotropic media and to the principle of maxi- 
mum energy release rate for smooth curvilinear cracks in anisotropic media. It can 
also be interpreted physically in this limit based on energetic considerations in the 
traditional framework of continuum fracture mechanics, in support of its general 
validity for real systems beyond the scope of phase-field models. Analytical predic- 
tions of crack paths in anisotropic media are validated by numerical simulations. 
Simulations also show that these predictions hold even if the phase-field dynamics is 
modified to make the failure process irreversible. In addition, the role of dissipative 
forces on the process zone scale as well as the extension of the results to motion of 
planar cracks under pure antiplane shear are discussed. 
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1 Introduction. 



The prediction of the path chosen by a crack as it propagates into a brittle 
material is a fundamental problem of fracture mechanics. It has classically 
been addressed in a theoretical framework where the equations of linear elas- 
ticity are solved with ze ro traction bou ndary conditions on crack surfaces that 
extend to a sharp tip (IBrobergl . Il999l ). In this description, t he stress distri- 
butions nea r the crack tip have the universal divergent forms ( iWillianj . 119571 ; 



Irwinl . I1957I ) 



(r,e) 



K 



/;?(0), 



where Km are the stress intensity factors (SIF) for the three standard modes 
I, II, or III of fracture (m = 1, 2 or 3), G is the angle between the radial vector 
of magnitude r with origin at the crack tip and the local crack direction and 
the explicit expressions of the /j/s are recalled in Appendix [C] (see Eq. flC.4|) . 
The allied energy release rate (or crack extension force) reads, for plane strain, 

G = a{Kl + Kl) + Kl/{2^i), (2) 



where v d enotes Poisson 's ratio, u is th e shear modulus, and a = (1 — z/)/ (2yu). 
Following Griffith f 192Cll ). Irwin ( 1957 ) postulated that for the crack to prop- 
agate, G must exceed some material dependent threshold Gc that is theoreti- 
cally equal to twice the surface energy (Gc = 27), but often larger in practice. 
Like other problems in fract ure, the prediction of the cra c k dire ction of prop- 
agation was first examined ( iBarenblatt and Cherepanovl . Il96ll ) for mode III 
which is simpler because the antiplane component of the displacement vector 
U3 is a purely scalar Laplacian field. In this case, the stress distribution near 
the tip, can be expanded as 



/i du3 K3 Q 
^30 = -^77 = /- — cos — - /iAa sm e . . . 
r dB ^/2^Tr 2 



(3) 



The dominant divergent contribution is always symmetrical about the crack di- 
rection. As a consequence, the knowledge of K3 al one cannot predict any other 
path t han a straight one. To avoid this impasse, IBarenblatt and Cherepanov 
(Il96ll ) retained the subdominant sin B term, which breaks this symmetry. 
They hypothesized that a curvilinear crack propagates along a direction where 
A2 = 0, when the stress distribution is symmetrical about the crack direction. 
In subsequent extensions of this work, several criteria have been proposed 
for plane loading, for which the tensorial nature of the stress fields makes 
it possible to predict non-trivial crack paths purely from the knowledge of 
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the s tress-intensity factors (iGoldstein and Salganikl . 1 19741 : ICotterell and Ricd . 
19801 ) ■ The generally-accepted condition = 0" assumes that the crack 



propagates in a pure opening mode with a s ymme trical stress distribution 

and is the direct ana- 
= for mode III. This 



19741 ) 



about its local axis (IGoldstein and Salganik 
log for plane strain (us = 0) of the condition A2 
"principle of local symmet ry" has been rationalized using plausible arguments 
flCotterell and Ricel . [T98oh but cannot be fully derived without an explicit de- 
scription of the process zone, where elastic strain energy is both dissipated 
and transformed nonlinearly into new fracture surfaces. As a result, how to 
extend this principle to anisotropic materials, where symme try considerations 
have no obvious generalization, is not clear (IMarderl . |200J). This is also the 
case for curved three-dimensional fractures although this appears little-noted 
in the literature. In addition, path prediction remains largely unexplored for 
mode III even for isotropic materials. 



Continuum models of brittle fracture that describe both short scale failure and 
macroscopic linear el asticity within a self-consistent set of equations have re- 
cently been propo s ed (lAranson et al.l. I2OOOI; 'Karm a et al.l . l200ll : lEastgate et al.l 



2OO2I : IWang et al.l . l2002l : iMarconi and Jagla . ,2005; ) . These models have already 
shown their usefulness i n vari ous numeric al simulations. For both ant iplane 
( iKarma and Lobkovskyl . |2004J ) and plane (IHenry and Levind . l2004j ) loading, 
they have proven capable to reproduce the onset of crack propagation at Grif- 
fith t hreshold as well a s dynamical branch i ng ins tabilities (IKarma and Lobkovskyl . 



2OO4J ) and oscillatory (IHenry and Levind . |200J) instabilities. In a quasistatic 



setting, this continuous media approach differs in spirit but has nonetheless 

much in common with a variational approach to brittle fracture (iFrancfort and Marigo , 



19981 ) proposed to overcome limitations of Griffith the ory. This is especiall y 



apparent w hen the latter is impleme r ited n umerically (iBourdin et al.l . I2OOOI ) 



using ideas ( Ambrosio and TortoreUi . 1990l ) initially developed for image seg- 
mentation (iMumford and Shahl . Il989l ). 



In this article, we analyze these self-consistent theories of brittle fracture 
and show how to derive laws of motion for the crack tip. This provides, in 
particular, relations which generalize the principle of local symmetry for an 
anisotropic material. Furthermore, we validate these relations by phase-field 
simulations. This validation is carried out both for the traditional variational 
formulation of the phase-field model with a so-called "gradient dynamics", 
which guarantees that the total energy of the system, i.e. the sum of the elas- 
tic and cohesive energies, decreases monotonously in time, and for a simple 
modification of this dynamics that makes the failure process irreversible. We 
find that both formulations yield essentially identical crack paths that are well 
predicted by the laws of crack motion derived from the phase-field model. 

For clarity of exposition, the relations derived from the phase-field model are 
summarized first in section |2] and interpreted physically in the context of 



3 



previous results from the fracture community. This section stresses why the 
second component of the Eshelby configurational force perpendicular to the 
crack axis is both physically meaningful and important for the determination 
of crack paths, even though this force has been largely ignored in the fracture 
mechanics literature since it was introduced. 

Our approach is applicable to a large class of diffuse interface descriptions of 
brittle fracture. However, for clarit y of exposition, we base our derivation on 



the phase-field model introduced by lKarma et al.l (I2OOII ). As recalled in section 



El in this description, the displacement field is coupled to a single scalar order 
parameter or "phase field" 0, which describes a smooth transition in space 
between unbroken (0 = 1) and broken states (0 = 0) of the material. We focus 
on quasi-static fracture in a macroscopically isotropic elastic medium with 
negligible inertial effects. Material anisotropy is simply included by making 
the surface energy '~f{6), dependent on the orientation 6 of the crack direction 
with respect to some underlying crystal axis. In section HJ we analyze the 
quasi-static motion of a crack, perturbatively for small departure from Griffith 
threshold {\G — Gc\ -C 1) and small anisotropy. The crack laws of motion are 
shown to be determined in a usual manner by solvability conditions, coming 
from translation invariance parallel and perpendicular to the crack tip axis. 

A different deri vation is provided in section [S] by generalizing Eshelby tensor 



(jEshelbyi . 119751 ) to phase-field theories. The particular case of motion under 
pure antiplane shear is then discussed in section [61 Our analytical predictions 
are compared with numerical phase-field simulations in section [7| where we also 
examine the sensitivity of the results to the irreversibility of the failure process. 
Our conclusions and some further perspectives of this work are t hen presented 



i n sec tion [HI Further information on the phase-field model of [Karma et al. 



(|2001[ ) is provided in Appendix [A| in the simple context of a stretched one- 
dimensional band. Details of some of our calculations are provided in the 
foll owing appendices [Bl and [Q A short version of this work has been published 



in ( Hakim and Karma . 2005f ) 



2 An overview of the physical picture and main results in the clas- 
sical fracture formalism 



In the formalism of continuum fracture mechanics, crack propagation has been 
traditionally analyzed by considering the crack extension force G defined by 
Eq. ([2]). This is a purely configurational force that points along the crack 
axis in the direction of propagation where GSl is the amount elastic energy 
released when the crack advances infinitesimally along this axis by a distance 
SI. When considering the propagation of a general curvilinear crack, however, 
it is necessary to consider the extension of a crack at some small infinitesimally 



4 



angle 66 with respect to its current depicted schematically in Fig. [H 

Physically, one would expect a configurational force, distinct from G, to be 
associated with the extra amount of elastic energy that is released if the crack 
propagates by 61 along this new direction, denoted here by t, as opposed to 
propagating the same distance along its current axis, denoted by xi. 



This additional force on the crack tip was considered by lEshelbyi (Il975l ). It can 
be interpreted physically as producing a torque on the crack tip that changes 
the crack propagation direction so as to maximize the elastic energy released. 
The force that produces this torque must act perpendicularly to the crack 
propagation direction and its magnitude is simply 

Ge = hm — (4) 

^ se^o 69 ^ ' 



where 6G is the difference between the crack extension force along the new 
direction and the old direction, i.e. along t a nd x\ in Fig. JH This torque is 



analogous to the well-known "Herring torque" (lHerringl . ll95ll . p. 143) acting on 
the junction of three crystal grains of different orientations in a polycrystalline 
material, with the main difference that Gq is a configurational force in the 
present fracture context while the Herring torque is produced by a physical 
force associated with the grain boundary energy, jgbiO), which is generally 
anisotropic. This force acts perpendicularly to each grain boundary segment 
at the junction of three grains with a magnitude d'^gb/dO. 



^2 




2 



Fig. 1. Schematic representation of an infinitesimal extension P1P2 of the crack 
of length 51 at and angle 59 measured with respect to the crack axis. The arrows 
pointing perpendicular to the crack denote the two analogs Gq and Gco of the 
Herring torque associated with the directional dependence of the crack extension 
force and the fracture energy around the crack axis, respectively. 

This analogy suggests that there should generally be two torques acting on 
the crack tip. The first, already mentioned, is Eshelby's configurational elastic 
torque Gg associated with the directional dependence of the crack extension 
force in reference to the local crack axis. The second is the physical torque 
associated with the directional dependence of the fracture energy, defined here 
by Gc{9), which should have a magnitude dGc{0)/dO = Gee by direct trans- 
lation of Herring's result for fracture. It follows that the balance of forces at 
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the crack tip should yield two conditions. The first is the standard condition 
of the classical fracture formalism associated with the balance of forces along 
the crack axis, G = Gc- The second is a new condition 



Ge = G,g = 27,, (5) 

which corresponds physically to the balance of the two aforementioned torques 
acting on the crack tip. While Ge pulls the crack in a direction that tends 
to maximize the release of elastic energy, Gee pulls the crack in a direction 
that minimizes the energy cost of creating new fracture surfaces. The second 
equality on the right-hand-side of Eq. only holds in some ideal brittle limit 
where the fracture energy is equal to twice the surface energy, defined here by 
7(6'), and 79 = d'-f{6)/ d6. We note that this ideal brittle limit is exact for the 
class of phase-field models analyzed in this paper but at best only approximate 
even for a strongly brittle material such as glass. The issue of the quantitative 
evaluation of Gee, however, should be kept separate from its role in crack path 
prediction that is our main focus in this paper. 

To see how this torque balance condition provides an explicit prediction for 
the crack path, it is useful to derive an expression for Ge by elementary means, 
directly from the definition of Eq. (HI) , instead of by evaluating an Eshelby-Rice 



type integral around the crack tip (iRicd . Il968l : lEshelbyl . Il975l ) , as done later in 



this paper (see section [5] and Appendix C); while both methods yield the same 
answer, the former is more physically transparent. For this purpose, we use 
the known expressions for the new stress intensity factors and K2 at the 
tip (corresponding to P2 in Fig. [1]) of an infinitesimally small k ink extension 



of length 61 of a semi-infinite crack (lAmestoy and Leblondl . Il992l ) . In the limit 



of vanishing kink angle, these expressions are given by 

Kl = Ki- ^K28e/2 + ... (6) 
k; = K2 + K,6e/2 + ... (7) 

to linear order in 66 independently of 61, where Ki and K2 are the stress 
intensity factors at the tip (corresponding to Pi in Fig. [1]) of the original 
straight crack. Using Eq. ([2]) with these new stress intensity factors to define 
G{66), we obtain at once that 6G = G{69) — G{0) = —2aKiK266, and hence 
using Eq. (jl]), that Ge = —2aKiK2. Substituting this expression for Ge in the 
torque balance condition we obtain the condition 

7^ ^ce le 



2aKi aKi ' 

where second equality only holds in the ideal brittle limit as before. In the 
isotropic limit where Gco vanishes, this condition reduces to the principle of 
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local symmetry which assumes that the crack propagates in a pure opening 
mode {K2 = 0). In contrast, for an anisotropic material, K2 is finite with a 
magnitude that depends both on Ki and the local crack propagation direction, 
i.e. Gee depends on the direction of the crack with respect to some fixed crystal 
axis in a crystalline material. For simplicity, we have restricted our derivation 
to a situation where linear elasticity is isotropic (e.g., hexagonal symmetry in 
two dimensions), but Eq. ([8]) could straightforwardly be extended to a more 
general situation where linear elasticity is also anisotropic. 



The recognition that the torque balance condition can be used to deter- 
mine the general path of a crack in a brittle material is the central result of 
this paper. This condition sheds light on the physical origin of the principle 
of local symmetry in the isotropic limit and shows how it can be generalized 
quantitatively to anisotropic materials. Although the co nfigurational fo rce per- 
pendicular to the crack tip was considered explicitly by lEshelbyl (119751 ). is has 
been largely ignored until recently. This is perhaps because the displacement 
of a small segment of crack perpendicular to itself, which one might naively 
expect to result from such a force, would appear unphysical and unrecon- 
cilable with the irreversibility of the fracture process. While such a motion 
is unphysical, it should be clear from the present considerations that all the 
torques acting on the crack tip, both the elastic configurational torque Gg and 
the physical torque Gc9 linked to fracture energy anisotropy, have been ob- 
tained solely from the consideration of an infinitesimal, physically admissible, 
extension of the crack at a small angle from its axis. In equating these two 
torques at the crack tip, the main assumption made is that the dynamics on 
the process zone scale is able to sample different possible microscopic states 
so as to permit local relaxation to mechanical equilibrium. 



There have been more recent attempts to incorporate the Eshelby elastic 
torque in the classical fracture formalism, where fracture surf aces are treated 



as mathemat i cally sharp boundar ies extending to the crack tip (lAdda-Bedia et al. 



I999I : lOleagal . I2OOII : iMardeii . l2004j ) . These theories, however, have not produced 



an explicit torque balance condition analogous to Eq. ([8]) that can be formally 
derived and tested. From this standpoint, the phase-field framework has the 
advantage of removing many of the ambiguities that arise when considering 
the motion of the crack tip in the classical fracture formalism. In the phase- 
field framework, a torque balance condition can be rigorously derived from the 
condition for the existence of a propagating crack solution that is spatially dif- 
fuse on the inner scale of the process zone, and must match smoothly to the 
standard solution of linear elasticity on the outer scale of the sample size. This 
condition reduces to Eq. (or Eq. for isotropic elasticity) in the limit of 
vanishing crack velocity, but contains additional contributions for finite crack 
velocity associated with dissipative forces on the process zone scale. 



Interestingly, the results of the phase-field analysis show that the component 
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of the dissipative force perpendicular to the crack tip vanishes for propagation 
in isotropic media because both the stress distribution and the phase field are 
symmetrical about the crack axis in this case. Consequently, within the phase- 
field framework, dissipative forces do not change the condition K2 = for 
crack propagation in isotropic media. For propagation in anisotropic media, in 
contrast, small velocity-dependent correction to the torque balance condition 
([8]) arise because this symmetry is broken. 



3 The KKL phase-field model of fracture 



Fracture is generally described in diffuse interface models (lAranson et al.l. 

2OOO I: iKarma et all I2OOII : lEasteate etaP . I2OO2I : IWang et all . l2002l : lMarconi and Jada 
2OO5I ) as a softening of the elastic moduli at large strains. This can be done 
purely in term of t he strain tensor but i t pro duces field equations with deriva- 
tive of high order (IMarconi and Jaglal . l2005l ). Here, we adopt the alternative 
approach of introducing a supplementary field 0, a scalar order parameter or 
"phase- field" , that describes the state of the material and smoothly interpo- 
lates between intact (0 = 1) and fully broken (0 = 0) sta tes. For defini t eness , 
we base our derivation on the specific model proposed by lKarma et al.l (120011 ) 
with energy density S, 



E = Spf{{dj(f)}) + g{(j)){Sstr •a 



(9) 



where dj = d/dxj denotes the partial derivative with respect to the cartesian 
coordinate Xj (j = 1,2,3) and £ strain is the elastic energy of the intact material. 
The equations of motion are derived variationally from the total energy of the 
system that is the spatial integral 

E= [ d^x £ (10) 



of the energy density. In the quasitatic case, these are 



5E 



d£ 

d[djU^] 
d£ 



_ d£_ 
d£ 



:i2) 



The three Euler-Lagrange Eq. (ITT!) for the cartesian components Uk of the 
displacement vector {k = 1, 2, 3) are simply the static equilibrium conditions 
that the sum of all forces on any material element vanish. The fourth Eq. (|T2i) 
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for is the standard Ginzburg-Landau form that governs the phase-field evo- 
lution, with X a kinetic coefficient that controls the rate of energy dissipation 
in the process zone, i.e. it follows from Eqs. fill I) and f[T^ that 



dE r f6E^ ^ 



Xld'x{ — ] <0 (13) 



In the simplest case of an isotropic elastic medium and isotropic 0, the phase- 
field and strain energy are simply. 



(14) 
(15) 



where Uij = [diUj + djUi)/2 is the usual strain tensor of linear elasticity. No 
asymmetry between dilation and compression is included since this is not nec- 
essary for our present purposes. The broken state of the material becomes 
energetically favored when £^strain exceeds the threshold £c and g{(j)) is a mono- 
tonically increasing function of that describes the softening of the elastic 
energy at large strain (^'(0) = 0) and produces the usual elastic behavior 
for the intact material {g{l) = 1,5''(1) = 0). In addition, the release of bulk 
stress by a crack requires the function g{(f)) to vanish faster than 0^ for small 
( t>, as recalled in Appendix El We therefore choose q((/)) = 40^ — 3 0^, as in 



( iKarma et al.l . l200ll : iKarma and Lobkovskyl . |2004J : iHenry and Levind . 120041 ) 



With these choices, the isotropic interface energy is equal to 



7o 



g{(j)) - 0.7165^2^^, 



(16) 



as sho wn in Appendix[X](Eq. (1A.22P ). by repeating the analysis of lKarma et al. 



(120011). 



In the present paper, we analyze the case of a phase-field energy £pf{{dj(l)}) 
without rotational symmetry which gives an anisotropic interface energy. A 
simple example used for concreteness and for the numerical simulations is 
provided by a simple two- fold anisotropy in the phase-field energy 

^vf = f (iV^r + e^^<P^2<p) (17) 



^ Note that with coordinates x', y' rotated by 7r/4 with respect to the x, y axes, the 
phase-field energy reads, £pf = ^(l + e/A){dx'(l))'^ + |(1 - e/4)(9y0)^, . 
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The interface energy of a straight fracture interface oriented at an angle 6 with 
the X-axis arises from the variation of the elastic and phase fields in a direction 
transverse to the fracture, namely with y) = (j)[—x sm{6)+y cos(^)]. There- 
fore, the only difference between Eq. (1141) and Eq. fll7l) in this one-dimensional 
calculation of the interface energy (Appendix is the replacement of n by 
k[1 — (e/2) sin 2^] in the anisotropic case. The allied anisotropic interface en- 
ergy thus follows directly from the isotropic expression (|T6l) and reads 

7(^) = 7o\/l-(e/2) sin 2^ (18) 

It reduces of course to the isotropic surface energy 70 of Eq. ( |T6l) in the e 
limit. 

With the specific energies of Eq. 0151) and fll7p . the variational phase- field 
equations read 

djWijg{^)] = 

+ ed,y(P] - g\<P){£strain " Q = -dt(p (19) 

X 

where aij is the usual stress tensor for an isotropic medium 

CTij = XukkSij + 2^Uij (20) 

Our aim in the following sections is to analyze the laws that govern the motion 
of a crack tip in this self-consistent description. 



4 Laws of crack motion as solvability conditions 

4^.1 The tip inner problem 

In the phase-field description, the obtention of laws of motion for a crack tip 
can be viewed as an "inner-outer" matching problem. The phase-field Eq. (1191) 

introduce an intrinsic process zone scale ^ = \J k/{2£c). The "inner" problem 
consists in the determination of a solution of Eq. f[T^ at the process zone scale 
^. The boundary conditions on this inner problem are imposed at a distance 
from the crack tip much greater than the process zone scale (r 3> ^ and much 
smaller than any macroscopic length. They should coincide with the short- 
distance asymptotics of the "outer" problem, namely the usual determination 
of the elastic field for the crack under consideration. Therefore, the imposed 
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boundary conditions on Eq. (fT9|) are that the material is intact (0-^1) away 
from the crack itself, and that for mixed mode I/II conditions, the asymptotic 
behavior of the displacement field is 

u^{r, 0) ~ ^) + ^2 4'(©; V)] (21) 



where /i is the shear modulus and the functions d™" are directly related to 
the universal divergent forms of the stress (Eq. [T]) and are explicitly given (in 
polar coordinates) by Eq. flC.6p and flC.7p of Appendix O The values of Ki 
and K2 are imposed by the boundary conditions at the macroscopic scale and 
do not significantly vary when the crack tip advances by a distance of order 
^. In the frame of the crack tip moving at velocity f , Eqs. ( [T9|) thus read 



9, [a,, (7(0)1 = (22) 

1) 

kV^0 - g'i(p){£ strain " £c) = d^cf) - eK dc,y(f) 



4-2 Perturbative formalism and solvability conditions 



Our first approach for obtaining the laws of crack tip motion consists in ana- 
lyzing the slowly moving solutions of Eq. fl22|) with boundary conditions fl2Tl) 
perturbatively around an immobile Griffith crack. For isotropic elastic and 
phase field energies and a pure opening mode, this Griffith crack corresponds 
to the stationary solution that exists for a{Kiy = Gc- Accordingly, we con- 
sider, for a small departure from Griffith threshold, 6K1 = \Ki — Kl\/ <^ 1 
and for a small K2 -C K^, a slowly moving crack with a small (two-fold) 
anisotropy in 0-energy (Eq. (ITSl) ). Our aim is to find the relations between K2 
and the anisotropy, as well as between Ki, K2 and the velocity required for 
the solution existence. 

Linearization of Eqs. fl22l) around the isotropic Griffith crack uf \ (f)^^'^ with the 
substitutions Ui = uf'^ + uf'\ = 0*^°^ + 0*^^^ gives. 



9,[a«^7(0W)] + d,[alfg'{<P^'^)<P^'^] = (23) 



11 



This can symbolically be written as 




X 



V 



\ 



I 











— e/t 











(24) 



where L is the linear operator on the left-hand-side (1. h. s. ) of Eq. fl23l) . The 
boundary conditions at infinity are that ^'-^^ vanishes and that u^^^ behaves 
asymptotically as in Eq. f l2Tl) but with K\ replaced by bK\^ the small departure 
from Griffith threshold, and is also assumed to be small. 

The linear operator L possesses two right zero-modes, that arise from the 
invariance of the zeroth-order problem under x and y translations, and can 
be explicitly obtained by infinitesimal translation of the immobile Griffith 
crack. For a general linear operator, the determination of the left zero-modes 
would nonetheless be a difficult problem. However, the variational character 
of the equations of motion imposes quite generally that L is self-adjoint (see 
Appendix [B]) and that left zero-modes are identical to right zero- modes. Thus, 
taking the scalar product of the two sides of Eq. fl2^ with the two translation 
zero modes provides two explicit solvability conditions for Eq. fl24l) . 

The scalar product with a left zero- mode (wf, U2,4'^) can generally be written 



Since the left vector is a zero mode of £, the only contribution to the 1. h. s. of 
Eq. (l25|l comes from boundary terms. 



fds n, {[^faj;) - «f VJ] ^(0(°)) + k^0« - u^U'] 9'i^^'^)^if + «:[0^5.0« - <P^'^d,<P^]} 



where n is the outward contour normal and the contour integral is taken 
counterclockwise along a circle (of radius r) centered on the fracture tip. 





(26) 
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4-3 Translations along x and crack velocity 



The zero mode corresponding to translations along x is {dxU^i \ dxU^2 \ dx4>''^^)- 
On the right-hand-side (r. h. s.) of Eq. fl25|) . the term proportional to the 
anisotropy e vanishes (by sjTiimetry or explicit integration). The r. h. s. of 
Eq. (l26l) can be simplified since on a circle of a large enough radius, g{(j)^^^) 
equals unity everywhere except in the region where the circle cuts the fracture 
lips. This region of non constant is far away from the crack tip where the 
crack is to a very good approximation invariant by translation along x and 
(9^.0 ^ ~ 0. Therefore, 



jj dxdy (dxuf\ dxU^ , dx 



Ki5Ki 

I 



l-v) 



(27) 



where the explicit formulas (]C.6|C.7p for the elastic displacements around 
a straight crack, have been used to obtain the last equality as detailed in 
Appendix O (see Eq. (1C.14I) ). Comparison between Eq. ( |271) and Eq. (1251) 
finally provides the natural result that the crack velocity is proportional to 
the departure from Griffith threshold. 



^jj dxdy[dx(, 



,(0)l2 



Ki5Ki 



'l-iy) = 6G 



(28) 



4-4 Translations along y and crack direction 



A second condition on crack motion arises from the zero mode correspond- 
ing to translations along y, {dyUi\ dyU2\ '^y4'^^^) ■ this case, only the term 
proportional to the anisotropy e contributes to the 1. h. s. of Eq. fl26|) . 

jj dxdydy<l)^'^^dxy(l)^"^^ = -j rfi/[9,0(°))U=_oo]' (29) 



Similarly to Eq. (l27j) . the r. h. s. of Eq. (126!) simplifies when the integration 
contour is a large enough circle 
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jj dxdy 



dylLl \ dyU2 \ dy 



: <l> ds riilu'f'^^a. 



(30) 



where again the exphcit evaluation in the last equality is detailed in Appendix 
[Cl (see Eq. f lC.20p . Thus, the second relation of crack motion reads 



+00 



^(l--) = f Jdy [dy<l>^'\-. 



(31) 



Eq. fl3T]) reduces to the principle of local symmetry (i. e. K2 = 0) for an 
isotropic medium and provides the appropriate generalization for the con- 
sidered anisotropy. Before further discussing these results and their physical 
consequences, we present a different derivation in the next section. 



5 Generalized Eshelby-Rice integrals 



The second approach, which we pursue here, directly exploits the variational 
character of the equations of motion and their invariance under translation. It 
yields identical solvability conditions as the approach of section H] when G — Gc 
and symmetry breaking perturbations are small, but it is more general since 
it does not require these quantities to be small. 



5.1 The Generalized Eshelhy tensor 



As shown by E. Noether in her classic work (INoetherl . Il918l ). to each con- 
tinuous symmetry of variational equations is associated a conserved quan- 
tity (charge) and an allied divergenceless current. Space (and time) trans- 
lation invariance are we ll-known to give the diverg e nless ene r gy-rn omentum 
tensor in field t heori e s (ILandau and Lifshita. 119751). lEshelbvl (119511) and fol- 



lowing authors (|Ricel. 119681: lEshelbvl. 119751 : iGurtin and Podio-Guiduglil . 119981 : 
Adda-Bedia et all . Il999l " lOleagal . 1200 ll ) have shown the usefulness of the anal- 



ogous tensor for classical elasticity theory. Here, we consider the generalized 
energy-momentum (GEM) tensor which extends Eshelby tensor for linear elas- 
tic fields (lEshelbvl . Il975h by incorporating short-scale physics through its ad- 
ditional dependence on the phase-field (p. 
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We find it convenient to define the four- dimensional vector field V'" = Ua for 
1 < a < 3 and ip"' = (p for a = 4, where the components of the standard 

displacement field. The inner problem Eq. (1221) can then be rewritten in the 
condensed form 

d£ d£ 

-5.,,vx-'d,<f> = d,^^^-^, a = l,...,4. (32) 



where here and in the following summation is implied on repeated indices 
(from 1 to 3 on roman indices and from 1 to 4 on greek ones). Chain rule 
differentiation provides the simple equality, 

= (33) 



Using Eq. (1321) to eliminate d£/dipa from the r.h.s. of Eq. (l33i) . we obtain 
djTij = -di(j)d,(f) for i = 1,2. (34) 



where the generalized energy-momentum (GEM) tensor T^j reads 



The GEM tensor T g is the sought extension of the classical Eshelby tensor 
Tj^ (lEshelbyl . Il95ll ) of classical linear elasticity 



Tj^j — £strain^ij (^jk^iUk (36) 



The GEM tensor Tij reduces identically to T^^ in the intact material where the 
phase-field is constant (0 = 1). Both tensors are non-symmetric in their two 
indices. The divergence of the GEM tensor taken on its second indice vanishes 
in the zero- velocity limit, when dissipation in the process zone also vanishes. 



5.2 Laws of crack motion 



In order to take advantage of Eq. (I34p . we integrate the divergence of the GEM 
tensor over a large disk Q centered on the crack tip (see Fig. [2]), following 
Eshelby computation of the conf igurational for ce on the crack tip treated as 
a defect in a linear elastic field (lEshelbyl . Il975l ) and subsequent attempts to 
derive criteria for crack propagation and stability (iGurtin and Podio-Guiduglil . 
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19981 : lAdda-Bedia et al.l . Il999l : lOleagal . l200ll ). The important difference with 
these previous computations is that, here, the GEM tensor (135!) is well-defined 
everywhere, so that the crack itself is included in the domain of integration. 
The integral of the divergence of the GEM tensor can be written as a contour 
integral over the large circle dQ bounding the disk Q, 



Fi = Jds Tij rij + jds Tij rij j dx di 

Ca^b b^a n 



0. 



(37) 



We have decomposed the circle dVL into: (i) a large loop Ca^b around the 
tip in the unbroken material, where A {B) is at a height h below (above) the 
crack axis that is much larger than the process zone size but much smaller 
than the radius R of the contour, ^ ^ /i ^ -R, and (ii) the segment {B — A) 
that traverses the crack from B to A behind the tip, as illustrated in Fig. O 
In both integrals, ds is the contour arclength element and rij the components 
of its outward normal. 

/ X2 £1 \ 





X2=+h 






X2=-h 



R 



Fig. 2. Spatially diffuse crack tip region with cp = 1/2 contour separating broken 
and unbroken material (thick solid line). 

Eq. (|371) provides an alternative basis to predict the crack speed and its path 
for quasi-static fracture. The Fj's can be interpreted as the parallel (i = 1) 
and perpendicular [i = 2) components with respect to the crack direction, of 
the sum of all forces acting on the crack tip. In Eq. fl371) . the three integrals 
terms from left to right respectively represents configurational, cohesive, and 
dissipative forces. We examine them in turn. 



5.2.1 Configurational forces and Eshelby torque 

We take A and B far back from the tip and close to the crack on a macroscopic 
scale but with the distance h between A and B much larger than the process 
zone scale. Namely, we consider the mathematical limit h +oo, R +oo 
with h/R ^ where R is the distance from A and B to the crack tip. In this 
limit, the first integral in Eq. (1371) is taken on a path that is entirely in the 
unbroken material where is constant and equal to unity. Thus, the tensor 
Tij reduces to the classical Eshelby tensor T^^ (Eq. (136|1 ) the first integral in 
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Eq. fl37|) yields the two components of the usual configurational forces p^'^°'^f\ 
Ft'"'^^= J dsT^^n, (38) 

Ca^b 



The first component, i^^^™"-^-*^ is the crack extension force and also Rice's J 



integral ( IRicd . Il968l ) . 



p(conf) ^ dsTf^n,, (39) 



Ca- 



With the known forms of the elastic displacement fields near the crack tip, as 
detailed in Appendix[C](see Eq. flC.281) ). one obtains the well-known expression 
([2]) of the crack extension force, 

Fi""""^^ = G = a{Kl + Kl), (40) 



The second component F2^°^'^^ can be computed in an analogous way from the 
elastic displacement fields near the crack tip, (Eq. (]C.29P and one obtains 

piconf) ^ _2aK^K2, (41) 



As discussed earlier, ^2^™"-^^ is the Eshelby torque ( Eshelby . 1975 ) that is 



readily interpreted physically if one imagine extending the crack tip by a 
small amount at a small angle 6 from the main tip axis. Then F2'^°^^^ is equal 
to the angular derivative of the crack extension force G{9) at = 0. 



This equality can be seen in two ways. First, we can use the general properties 
of Eshelby tensor. We denote with a tilde the elastic quantities corresponding 
to the crack with the small extension of length s at an angle 9. Since the crack 
extension is along the direction (cos(^), sin(^)), we consider the allied vector 
obtained from the Eshelby tensor, T^^ = cos{9)Tfj + sm{9)T2j. The flux of 
this vector vanishes when taken through the contour that goes along the great 
circle from A to B and then continues in a classical way along the lips of the 
extended crack, as drawn in Fig. [3] 

+ /+/+/ + /)rfsf,>, = (42) 

Ca^b B~^C C^D Cd-.e E^F F^A 



The two integrals on the fracture lips from C to D and E to F do not con- 
tribute since the integrand vanishes: the 9 direction is along the path and the 
normal stresses vanish on the fracture lips. The same argument shows that 
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Fig. 3. Sketch showing the contour integral decomposition in Eq. (j42p . The crack 
with its virtual extension at an angle 9 is depicted by the the thick bold line. The 
integral contour follows the great circle from A to B; it continues along the upper 
lip of the crack from B to C and then along the upper lip of the virtual extension 
from C to D; it then encircles the extended crack tip following the small circle from 
D to E; finally it comes back to A along the lower crack lips via E and F. 

the integrand is simply equal to ±Sstrain sin(6') for the integrals from B to C 
and F to E along the lips of the original fracture. The integral on the small 
circle around the extended crack tip Cd^e is equal to —G{6) where G{6) is 
the energy release rate at the end of the small crack extension. Eq. ( H2l) thus 
reduces to 



j [cos{e)f^^ + sHe)fi^ n, = G{e) + J rfxsin(^^)[4t»„(x) - 4™.n(a;)](43) 

Ca^b R 

where £ttrain ^iid S^train respectively denote the elastic strain energy densities 
on the upper and lower fracture lips. The required identity between F2°^^'^ and 
the angular derivative of dG/ d9\g=Q follows from differentiation of Eq. ( l43l) with 
respect to 6' at ^ = 0. When this is performed, there are two kinds of terms. 
Terms coming from the differentiation of the explicit trigonometric functions 
in Eq. (H3l) and terms coming from the implicit dependence upon 9 of tilde 
quantities. However, in the integral on the l.h.s of Eq. (H3l) . these implicit terms 
vanish as the length s of the extension is taken to zero, and in the integral on 
the r.h.s. they are multiplied by a vanishing sine function. Moreover, for the 
straight fracture at 6* = 0, only Gxx is non-zero on the fracture lips and it is 
of opposite sign on the upper and lower fracture lips (Eq. flC.4l) ). The strain 
energy densities which are quadratic in the stress a^x are therefore equal on 
the upper and lower fracture lips and after differentiation, the contribution 
of integral term on the r.h.s of Eq. fH5]) vanishes at zero [fj. Finally, in the 

^ Subdominant terms, coming for instance from a macroscopic curvature of the 
crack, could be different on the two crack lips but note that the integral range is on 
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limit of a vanishing extension length (s — 0) tilde quantity tend toward their 
(non-tilde) values on the original fracture and one obtains 

^(™"/) _ J rrE^ _ dG{e) ^ 

Ca->b 



Fr'>= I T2>, = lhn^|,=o = G,(0) (44) 



This relation between the second component F2^°'^'^^ and the angular derivative 
of G{6) can also be obtained by comparing their explicit expressions in term 
of the SIF Ki and K2. As it is well known, the SIF Ki and K2 at the end 
of a small extension can be expressed as linear combination of Ki and K2 



(lAmestoy and LebloncO . Il992l ) 



k^ = F^,i9)K^ + F^2{9)K2 

k2 = F2i{e)K^ + F22{e)K2 (45) 

with clearly Fii(O) = ^22(0) = 1 and Fi2(0) = F2i(0) = and the derivative 
at 6* = 0, F[^{0) = FU9) = 0, as alreadv mentionned in section[2]( Eq. flHllTD ). 



A detailed computation ( lAmestoy and Leblondl . Il992l ) provides the other two 



derivatives F[2{Q) = —3/2 and -^21(0) = V^- Therefore, one obtains 

lim = 2aKiK2[F[2{0) + F^^{0)] = -2aK^K2 (46) 

s^o dU \e=o 



This is indeed identical to the expression of F2^°^^'^ obtained by a direct com- 
putation (Eq f HT]) ) and it provides a second derivation of Eq. (jH 



5.2.2 Cohesive forces 

An important new ingredient in Eq. ( 1371) is the second portion of the line 
integral (/g^^) of the GEM tensor that traverses the crack. This integral 
represents physically the contributions of cohesive forces inside the process 
zone. To see this, we first note that the profiles of the phase-field and the 
three components of the displacement can be made to depend only on X2 
provided that the contour is chosen much larger than the process zone size 
and to traverse the crack perpendicularly from B to A. With this choice, we 
have that rii = —1, n2 = 0, along this contour and therefore that, for i = 1 

+h 

Fi^"^^ = J ds Ty nj = - J dx2Tu (47) 



a length scale that is vanishingly small on a macroscopic scale 
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The spatial gradients parallel to the crack direction (diip^) give vanishingly 
small contributions in the limit h/^ ^ +00 and R/^ +00 with h/R ^ 0. 
Thus, the integrand on the r.h.s of Eq. fl47p reduces to the energy of a Id 
crack which, as recalled in Appendix |X] (Eq. (1A.22P ) is itself independent of 
the strain and can be identified to twice the interface energy 7 

+h 

= - J dx2S{<j), ^20, d2U2) = -27 (48) 

-h 



This yields the expected result that cohesive forces along the crack direction 
exert a force opposite to the crack extension force with a magnitude equal to 
twice the surface energy. 

One similarly obtains for i = 2, in the same limit ^ <^ /i <^ -R, the other 
component F2^°^'' of the force perpendicular to the crack direction 

+h +h +h 

Ft"^ = I dsT^.n, = - I dx2T2i = I dx2^^d,^P^ = J dx^^d^^Pm 

B~*A -h -h -h '''^ 



The last equality comes from the fact that the only considered anisotropy is 
in the phase field part Epf of the energy density and that, as above, gradients 
parallel to the crack direction give negligible contributions far behind the crack 
tip. F2^"^^ can be expressed as the angular derivative of the surface energy at 
the crack tip direction 6 = 0. For a material broken along a line lying at a 
direction 6 with the x-axis, the displacement and phase fields only depend 
on the normal coordinate rj = —xi sin(^) +X2 cos{6). The local energy density 
S[(j), di(j), (920, dr^Urj] is therefore equal to S[(j), — sin(6')9^0, cos(6')(9^0, drjUr^]. The 
allied surface energy reads 

+00 

2-f{9)= [ dri £[-sm{e)dr^^,cos{9)dr,(l),dnUn] (50) 



Differentiation with respect to 6 brings on the r.h.s. of Eq. (!50|) terms coming 
from the explicit dependence of the integrand on 6 as well as terms com- 
ing from the implicit dependence of the fields on the breaking angle (for an 
anisotropic material). However, the contribution of the implicit terms vanishes 
since for any given angle the fields minimize the total energy and no field vari- 
ation leads to an energy change at linear order. Therefore, one obtains 

—00 —00 
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since rj reduces to 0:2 for 6^ = 0. Comparison of Eq. (H9|) and (l5Ti) shows that 
F^'^ = -2^7|.=o (52) 



as announced. 

Of course, the relation fl52l) can also be checked by direct computation for 
any explicit form of the phase field energy. For instance, in the simple case of 
Eq. ( fT7|l . one obtains from Eq. ( I49l) 

F^"'^ = J dx,^d,<P =y dxMd2<Pr = £70 (53) 



where, for the last equality, it should be noted (see Appendix that the 
second integral in Eq. flS51) is equal to the energy (by unit length) of the cracked 
material which is itself equal to 270. The result of Eq. fl5^ indeed agrees with 
Eq. (152!) . e7o = —27^(0), since the interface energy in the direction 9 is given 
by Eq. ([HD. 

The force F^; '^°^^ is the direc t analog of the Herring torque 70 = ^7/ d6 on grain 



boundaries (IHerringl . Il95ll . p. 143). This torque tends to turn the crack into 



a direction that minimizes the surface energy. 
5. 2. 3 Dissipative forces 

The last term in Eq. (I57|) gives the two components of the dissipative force 

+00+00 

pidi^) ^ ^^-1 f f dxidx2 di(f)di(j), (54) 



The limit where the disk area f2 tends to infinity has been taken since the 
integrand vanishes outside the process zone. In contrast to the configurational 
and cohesive forces, the dissipative force clearly depends on the detail of the 
underlying diffuse interface model. 



5.2.4 Force balance and anisotropic generalization of the principle of local 
symmetry 

Substituting the results of Eqs. (139|) to (15^ into Eq. (!37j) . the two conditions 
of Eq. (1371) can be rewritten in the compact form 
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Fi = G - a - f}''''^ = 0, (55) 

F2 = G,(0)-a^^(0)-Ff^) = 0, (56) 

where we have used the fact that Gee = '^le- Eq. ( l55l) together with Eq. ( 15^ 
predicts the crack speed for G close to Gc 

^ {G-Ge) (57) 



jj dxidx2{di(f)o 



where 0o is the phase- field profile for a stationary crack (IKarma and Lobkovskyl . 



2004J ). and thus the integral in the denominator above is just a constant of 
order unity. Eq. fl56l) . in turn, predicts the crack path by imposing K2 at the 
crack tip, 

K2 = - (GeeiO) + Ff /{2aK,). (58) 



The component F2 *^ of the dissipative force vanishes with the crack velocity 
in the quasitatic limit. So, in this limit, the microscopic details of the process 
zone do not play a role and the crack direction is uniquely determined by the 
directional anisotropy of the material through the simplified condition 

K2 = -Gee{0)/{2aK^). (59) 



Eq. fl59l) replaces the principle of local symmetry for a material with an 
anisotropic surface tension energy. It reduces of course to the principle of 
local symmetry in an isotropic material, since then Gee vanishes. One can also 
note that quite remarkably, Eq. (159|1 only contains macroscopically defined 
parameters and is independent of the detailed physics of the process zone. 

Outside the quasistatic limit, K2 = should continue to hold for an isotropic 
material since ^2^'^*'^'' vanishes even for a finite crack speed. The latter follows 
from the symmetry of the inner phase-field solution for a propagating crack 
with K2 = 0, 0(xi,a:2) = (f>{xi, —X2), which implies that the product di(j)d2(p 
in Eq. (|54|) is anti-symmetric and that the spatial integral of this product 
vanishes. In an anisotropic material, however, (p is generally not symmetri- 
cal about the local crack axis and ^2^'^*'^^ should generally be non-zero. The 
crack direction should then become dependent on the details of the energy 
dissipation in the process zone. 

A small velocity expression for the dissipative force perpendicular to the crack 
axis can be obtained by considering the phase-field profile that corresponds 
to a stationary Griffith crack in an anisotropic material. Here, (pQ, is uniquely 
defined as the stationary phase-field profile that exists for a unique pair of 
values of Ki and K2 that satisfy the conditions of equilibrium parallel and 
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perpendicular to the crack axis, a{Kl + K2) = Gc{0) and —2aKiK2 = Gc9{0), 
respectively. For small velocity, Eq. fl5^ must therefore reduce to Fg^'^*''-* = 
fX~^/(0) where the integral 

/(O) = [ dXidX2 ^10^^20^ (60) 



is a dimensionless constant that, like Gc and Gee, depends generally on the 
local orientation of the crack with respect to some fixed reference axis chosen 
here as ^ = 0. For small velocity, Eq. fl58|) therefore becomes 

K2 = - {G,e{0) + vx-'l{0) /(2ai^i), (61) 



where /(O) vanishes in the isotropic limit since 0^ approaches 0o and hence 
becomes symmetrical about the crack axis in that limit. 



5.3 Comparison with the maximum energy release rate criterion 



The principle of local symmetry and the maximum energy release rate criterion 
gives slightly different results in general, for instance for the prediction of 
the finite angle of a kink extension at the tip of a crack. The two criteria 
coincide however for smooth cracks. It is interesting to note that it is also true 
for the present anisotropic generalisation (Eq. dSHD) of the principle of local 
symmetry. One way to generalize the maximum energy release rate criterion 
for anisotropic material is to require t he crack growt h to take place in the 



direction that maximizes G{6) — 2'y{6) (ILeblondl . l2005l ) where as before G{9) 
is the energy release rate for an infinitesimal extension at the crack tip at 
an angle 6 (where as before ^ = is the direction of the unextended crack). 
For a smooth crack, the condition that this quantity be maximal in the crack 
direction yields 

±[G{e) - 2^{e)]e=o = (62) 



With the help of Eq. fH6|) . this is seen to be identical to Eq. fl59l) as stated. 



5.4 Crystalline materials 



Our results have interesting implications for crack propagation in crystalline 
materials. Basic experimental studies have demonstrated the existence of both 
"cleavage cracks", which are cracks that propagate along low energy crystal 
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silicon, and smooth cracks (jPeegan et al. 



planes, such as {111} (iHauch et all Il999 ) or 1 1 10) ( jPeegan et al.l . l2003l ) in 



2003h that resemble qualitatively 



the cracks seen in isotropic materials. While the propensity for crack propaga- 
tion along cleavage planes in crystalline materials is to be expected energeti- 
cally, the observation of smooth cracks in those same materials is perhaps less 
intuitive. Theoretical attempts have been made to understand when cracks 
will cleave crystals using both en ergetic arguments and lattice simulations 



( jPeegan et al.l . l2003l : iMarderl . 120041 ). However, a consistent theoretical picture 



has not yet emerged. 

The crack propagation law Eq. (1591) provides an explicit prediction of when a 
crack will propagate along a cleavage plane, or smoothly in other directions. 
Restricting our discussion to two dimensions for simplicity, the surface energy 
in a crystalline material is expected to show a cusp behavior 



7(0) = 7o(l + 



+ 



(63) 



near a cleavage plane (and more generally near sets of equivalent low energy 
crystal planes imposed by symmetry), where 9 measures the angle of the sur- 
face away from this plane, and to be a smooth differentiable function of 9 for 
other orientations. In terms of the physical picture outlined in Section 2, this 
cusp behavior implies the presence of a finite Herring torque on any small 
extension of a crack at an infinitesimal angle away from a cleavage plane. 
Therefore, a crack will be essentially trapped along a cleavage plane until the 
Eshelby configurational torque is large enough to tilt the crack away from this 
plane. Restated in terms of the propagation law, Eq. fl59l) can be obeyed for 
small non-zero angles only when \K2\ exceeds a threshold K'f' with 



for G^Gr 



(64) 



(c) 

Consequently, |i^2| should exceed K2 for a cleavage crack to change direc- 
tion. Eq. (159|) also implies that a crack will propagate smoothly for other 
orientations away from cleavage planes where the surface 7-plot is smooth. 

One interesting prospect to test this prediction is to examine its consequences 
for thermal fracture in crystalline materials, where quasistatic oscillatory cracks 
have been studied under well-controlled experimental conditions. In particu- 
lar, experiments have revealed that the onset of crack oscillations is markedly 
different in anisotropic and isotropic materials. In crystalline silicon wafers 
that cleave preferentially {110} planes, the onset of crack oscillations is de- 
layed in comparison to an isotropic material and is accompanied by a discon- 
tinuo us jump in oscillati on amplitude consistent with a subcritical bifurca- 
tion (jPeegan et al.l . 120031 ). In contrast, the onset of thermal crack oscillations 
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in isotropic material has been shown to be supercritical in a recent phase- 
field mo deling study, con sistent with earlier experimental observations in glass 
(see (j Corson et al.l . |2008| ) and earlier references to the experimental literature 
therein). The existence of a finite threshold Eq. fl64p to escape a cleavage 
crack precludes a smooth transition to crack oscillations around a cleavage 
plane. One would therefore expect a subcritical bifurcation for the onset of 
crack oscillations if Eq. fl59|l is used in conjunction with a typical 7-plot for a 
crystalline material that exhibits cusps. However, a detailed study is clearly 
needed to validate this expectation and to make contact quantitatively with 
experiments. 



6 Motion under pure antiplane shear 



As recalle d in the introduction, the pr i nciple of local symmetry was first pro- 
posed in (IBarenblatt and Cherepanovl . Il96ll ) for crack motion under pure an- 
tiplane shear. This particular case does not seem to have attracted much 
interest subsequently, presumably because rotation of the crack front is ob- 
served and fractures under mixed mode I-III load i ng are found to be unstable 
in three dimensional isotropic materials (jSommeii . Il969l ) . The criterion of mo- 
tion under antiplane shear could nonetheless have some importance for the 
development of the tridimensional instability. It is also conceivable that 2D 
motion under pure mode III loading could be effectively realized in an ap- 
propriate anisotropic material, like for instance a thin layer of sintered glass 
beads. We therefore find it interesting to briefiy examine this criterion for mo- 
tion under pure antiplane shear with the formalism developed in the previous 
sections. 



Since in a pure mode III motion, the displacement field reduces to its third 
component U3 that is a purely scalar Laplacian field, the diverging stress dis- 
tribution near the tip Eq. ([3]) is always symmetrical and produces no configu- 
rational force perpendicular to the crack axis. Consequently, for propagation 
in an anisotropic material, the propagation law reduces simply to the condi- 
tion that the Herring torque vanishes, 75 = 0. This condition implies that in 
the limit of vanishing velocity, a quasistatic crack propagates in a direction 
that corresponds to a local minimum of the surface energy; it can be argued 
that propagation in a direction of maximal 7 is unstable because the config- 
urational torque amplifies small departures from this direction. Furthermore, 
for finite velocity, the dissipative force perpendicular to the crack axis, F2'^^\ 
also vanishes since the phase-field profile must be symmetrical about the crack 
axis for a direction where 75 = 0. 

For propagation in an isotropic material, the situation is more subtle than for 
the mode I/II case. The evaluation of the contour integral that is the direct 
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analog for mode III of the r.h.s. of Eq. (130|1 (or equivalently Eq. (l38ll ). gives 
only a no n- vanishing force perpendicular to the crack axis if the sub dominant 
antisymmetrical contribution of the stress distribution (i.e., the second term 
on the r.h.s. of Eq. is included. This force is proportional to K3A2VR, 
where R is the radius of the integration contour around the crack tip, where 
the square-root behavior follows from dimensional analysis. This force vanishes 
if ^2 = 0, thereby suggesting that the original formulation of the principle of 
local symmetry for mode III might be applicable. However, our inner-outer 
matching procedure used to compute this force is predicated on choosing R 
much larger than the scale of the process zone but vanishingly small on the 
outer scale of the system size set by material boundaries. Therefore, the mag- 
nitude of this force is left undetermined in the present analysis. Further work 
is therefore needed to determine if the inner and outer scales can be clearly 
separated for pure antiplane shear and if ^2 = can rigorously serve as a local 
condition to predict crack paths in isotropic material. 



Additional insight into this question can be gained by repeating the analysis 
of Section 2 for a small extension 61 of a mode III crack. The analogous 
expression for the stress intensity factor at the tip o f the extended crack is 
= — bfiA2\/6l66 to linear order in 69 f SihI . Il965l ) where 6 is a numerical 
constant, and hence Ge{0) ~ K3A2V6I. One important difference with plane 
loading is the square-root dependence of Ge{0) on the crack extension length, 
which is also reflected in the a/R dependence of the integral just mentioned 
above, which yields the configurational force perpendicular to the crack tip 
for mode III. Since the only natural cut off for the crack extension length on 
the outer scale of the system is the size ^ of the process zone, this result seems 
to imply that Gei^O) ~ K3A2\/^ up to a numerical prefactor. It also yields the 
local symmetry condition A2 = in the isotropic limit, where the symmetry 
of the phase-field profile makes the dissipative force Fa"^*^-* vanish. 



7 Numerical simulations and tests 



We focus here on testing the relation (13^ between K2 at a crack tip and the 
derivative of the interface energy, in the case of plane strain. We numerically 
compute the extension of a preexisting straight crack as described by the 
phase-field equation f|T9|) and f|T8|) . For a pure mode I loading of the preexisting 
crack and an anisotropic surface energy, Eq. fl59l) predicts that the growth of 
a kinked extension takes place at an angle 6 such that K2 is adequate at the 
growing tip. More explicitly, on the one hand, Eq. fHS]) gives K2 at kink tip as 

K2 = F2i{e)KiC:^K,^ (65) 
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where the second equahty is vahd for small angles. Therefore, one obtains for 
small angles, 

9 

- 2akik2 ^ aK^ (66) 



since Ki is equal to Ki at dominant order in 6 (Eq. ( H5i) . On the other hand, 
the surface energy (ITSl) gives 

- 27e(0) = £70 (67) 



Eq. ( l59l) . which translates in the equality of the l.h.s. of Eq.( l66i) and ( 1671) . 
simplifies for G close to Griffith threshold when aKf ~ 270. Then, it simply 
gives for the initial angle 6 of the kink crack 

= { (68) 



which is strictly valid for e ^ 1 in the limit G Gc 



This prediction was tested numerically. Eq. (1191) was solved by using an Euler 
explicit scheme to integrate the phase-field evolution and a successive over re- 
laxation (SOR) method to calculate the quasi-static displacement fields Ui and 
U2 at each time step. We used as initial condition a straight horizontal crack 
of length 2W centered in a strip of length AW horizontally and 2W vertically, 
with fixed values of Ui and U2 on the strip boundaries that correspond to the 
singular stress fields defined by Eq. ([1]) for prescribed values of Ki and 7^2- 
We used A//i = 1 [a = 3/(8/i)], Sc/ fi = 1/2, a grid spacing Axi = Ax2 = 0.1^, 



and W = 50^, where the process zone size ^ = ^Jk/{2Sc)- We checked that 
the results were independent of width and grid spacing. 

We first verified that, in the isotropic limit, the kink angle was well predicted 
by the local symmetry condition K2 = , which implies that 6 ~ —2K2/K1. 
Then, for the anisotropic case, we chose K2 = and G just slightly above Gc- 
The results for the kink angles observed for several simulations with different 
magnitudes of the surface energy anisotropy e are shown in Fig. |H The pre- 
diction (|68ll is seen to be in good quantitative agreement with the results of 
the phase-field simulations. 

We have also tested our prediction for pure mode III cracks. We used the 
same phase-field model and anisotropy form of 7(6') as for the plane strain 
case, albeit with the strain energy corresponding to pure antiplane shear 



^strain 2 



{diu.f + {d2U,f . (69) 
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Fig. 4. Kink angle 9 versus surface energy anisotropy e for plane strain predicted as 
6 = e/2 and extracted from phase- field simulations (filled circles) for G/Gc « 1.1. 
Inset: phase- field simulation showing = 1/2 contours equally spaced in time for 
e = 1.2. 

We used as initial condition a straight horizontal crack of length ?)W centered 
in a strip of length QW horizontally (along the Xi axis) and 2W vertically 
together with Ed = 1/2, the anisotropy e = 1.8, a grid spacing Axi = Ax2 = 
^/6, a half strip width W = 50^, and a fixed displacement = 11.313 ^ on 
the X2 = i^W boundaries corresponding to a crack slightly above the Griffith 
threshold, where as before the process zone size ^ = ^k/ {2Sc)- 

The results of this simulation shown in Fig. [S] confirm that a crack centered 
initially in the strip with its axis parallel to the ^ = direction, kinks at a 
45° angle {6 = vr/4) that is consistent with the analytical prediction 70 = 
for mode III in an anisotropic material. 

Fracture in the phase-field model that we have considered here is a reversible 
process in the sense that cracks can (and do) heal when stresses are removed. 
This can also be observed in some experiments under very clean conditions 
when no alterations of exposed surfaces follow breaking. Nonetheless, this is 
sometimes considered a troublesome feature since it does not occur in more 
usual conditions. One could think of introducing irreversibility in a "physical 
way" by adding another field to mimick surface oxydation. In a simpler but 
more ad-hoc fashion, one can only accept evolutions that decrease the value 
of the phase field. To assess the importance of reversibility on our results, the 
numerical simulations above for mode I /II and mode III were redone with this 
second scheme ((i.e. taking Eq. (fT9|) as written when dt4> < and otherwise 
replacing it by dt(j) = 0). Reassuringly, the results are essentially identical to 
those plotted in Fig. H] where the 0=1/2 phase-field contours superimpose 
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Fig. 5. Phase-field simulation for pure antiplane shear for e = 1.8 (</> = 1/2 contours 
are equally spaced in time). The dashed line is a guide to the eye for the 45° kink 
angle predicted by the vanishing torque condition 70 = 0. 

perfectly for the two sets of simulations with and without reversible dynamics. 
This insensitivity of the results to the introduction of irreversibility does not 
appear surprising since our derivation of crack propagation laws for modes 
I/II and III rely on the existence of propagating solutions for which dtcf) < 0. 



8 Conclusion 

We have analyzed here the laws of quasistatic crack tip motion within the 
phase-field framework. The analysis provides a derivation of the principle of 
local symmetry and of its generalization to anisotropic materials. It also un- 
derlines the role of the configurational force perpendicular to the crack tip 
direction. The results can be interpreted physically as a simple force balance 
condition. The variational character of the phase-field equations of motion 
played an important role in the derivation of the equations for the crack tip. 
It directly allowed us in section [5] to define a generalization of Eshelby tensor 
that includes the phase-field and short-scale physics while keeping its diver- 
genceless property. Its role in the derivation of section H] may be less central 
from a conceptual point of view but the self-adjointness of the linear operator 
allowed the obtention of explicit formulas. In any case, it should be noted that 
direct link to energy considerations a la Griffith require a variational model. 

Several questions appear worth of further investigations. First, even if pure 
mode III cracks are not realized under most experimental conditions, the appli- 
cability of the principle of local symmetry {A2 = 0) remains to be established 
on a firmer footing in the isotropic limit where the force perpendicular to the 
crack tip seems to depend on an arbitrary cut-off on the scale of the process 
zone. Second, the crack propagation laws have been derived for a gradient for- 
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mulation of the phase-field dynamics where failure is reversible. Even though 
our numerical simulations indicate that crack paths are not altered by a simple 
ad-hoc introduction of irreversibility where the phase-field can only decrease, 
the role of irreversibility is worth investigating more fully. 

The phase field energy is certainly another aspect that would benefit from 
further refinement. As written in Eq. fll5lll9p . it does not distinguish be- 
tween compressive and extensive strains which is quite an unphysical feature. 
Note tha t this property is also shared by the variational formulation advo- 
cate d by Francfort and Marigol ( 19981 ) (illustrated in numerical simulations 



of (IBourdin et al.l . l200d . Fig. 4) in which it is ref erred to as sample inter- 



prenetration). Some remedies have been proposed (IHenry and Levind . |200J) 
that break the variational character of the equation of motion and do not ap- 
pear entirely satisfactory. The development of more physically motivated and 
material adapted energies certainly appear as an interesting future endeavor. 

The extension of the present analysis to three dimensions where fracture paths 
are geometrically more complex is another important future direction. Nu- 
merical simulations and preliminary analysi s addressing this question will be 



reported elsewhere (jPons and Karmal . l2008l ). 



Finally, we hope that the results reported here will contribute to stimulate 
further experimental investigations of quasistatic crack motion of cracks in 
anisotropic media. 

We thank M. Adda-Bedia, J. B. Leblond, A. Chambolle, G. Francfort and J. J. 
Marigo for valuable discussions and instructive comments. A.K. acknowledges 
the support of DOE Grant No. DE-FG02-07ER46400 and the hospitality of 
the Ecole Normale Superieure in Paris where part of this work was completed. 



A The KKL phase-field model in one dimension 



In this appendix, we recall the analysis of the KKL phase-field model in one 
dimension, that is t he snap-back of a stretched elastic band, as described in 
( iKarma et al.l . l200ll ). In particular, the energy of the fractured state Eq. flA.221) 
provides the expression of the interface energy given in the main text. We also 
show how the fractured solution appears in one dimension in this model. For 
an elastic band of size 2L, the elastically stretched state is the only allowed 
state when the total strain 2A is low enough. Above a critical total strain 2Ac 
two other non trivial solutions appear via a saddle-node bifurcation, one be- 
ing dynamically stable and the other being unstable. At the bifurcation, both 
solutions have a higher energy than the elastically streched state. However, 
the dynamically stable solution becomes energetically favored as compared 
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to the elastically streched state when the total strain becomes higher than 
2Ag (> 2 Ac), which corresponds to Griffith threshold in the model. This sce- 
nario is illustrated by numerical solution in Fig. lA.li The unstable solution 
corresponds to the energy barrier (the Eyring state) that has to be overcome to 
create the fractured state and it provides the corresponding activation energy. 

For a one-dimensional band, the KKL energy reads. 




(A.l) 



with the function g monotonically increasing from g{0) = in the fully broken 
state to g{l) = 1 in the intact material with also g'{l) = to recover linear 
elasticity. Steady state solutions obey the equilibrium equation obtained by 
variation of Eq. flA.ip . 



dy[g{^)dyu] = 



(A.2) 
(A.3) 



with the boundary conditions u{±L) = ±A, (f){±L) = 1. The elastically 
stretched band 0=1, u{y) = yA/L is always a solution of Eqs. ( 1A.21 IA.3I) 
and its energy is equal to the usual purely elastic one 



E = (A + 2/i)AVL. 



(A.4) 



In order to analyze the existence of other less obvious solutions of Eqs. (lA.2IIA.3p . 
it is useful to note that Eq. (1A.3P can be integrated once to obtain 



(A.5) 



with c a constant yet to be determined. Eq. (1A.5I) allows the elimination of 
the strain field from the phase-field equation flA.2|) which then reads 



Kdyy(j) = Scg'{(j)) 



(A.6) 



In a usual way, it is helpful to consider y as a fictitious time and to think of 
Eq. (]A.6[) as describing the motion of a point particle in the effective potential 



VefA^) = ^^+g{^) (A.7) 
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With this analogy, a non-trivial solution of Eq. (1A.6P corresponds to a particle 
that starts at "time" y = —L from 0=1 with a negative velocity dy(j) < 0, to 
reach a minimum = 0^ at ?/ = where the "velocity" dy(f) vanishes; from 
this turning point it then follows the time-reversed motion and comes back to 
= 1 at ?/ = +L. The integrability of this one-dimensional motion gives the 
conservation law 

^(9,0)2 + K//(0) = l4//(0™) (A.8) 



This allows us to express the energy (Eq. ( lA.ll) of the corresponding non-trivial 
solution as 



E 



^"^^ I Ur .r ^J ^ + ^^fS^'^^) - '^am (A.9) 

L \J^eff{(t)m) - Veff[(t)) 



Two constraints determine the two unknown constants c and 0^ as a function 
of the dimensionless strip width I and dimensionless total strain 5. First, the 
particle motion should take a total time 2L with 

i = L^28Jk = I ci0 ^^^^^^ 

L V'^e//(0m) - 14//(0) 



Second, the overall integrated strain should equal the imposed total strain 
5 = A7(A + 2/i)/K = c / '^^ = (A.ll) 



We analyze more specifically the case of a macroscopic strip of width £ ^ 1. 

We find it convenient to first consider the dependence of Eq. (lA.lOp on c. Since 
the sum of kinetic energy [n/ {2£c){dy(t)Y)] and potential energy [Ve//(0)] is 
constant along the particle trajectory (Eq. (1A.8I) . the initial potential energy 
is always lower than the potential energy at the return point where the kinetic 
energy vanishes), \4//(l) < K//(0m)- For a given 0m, the time L spent by the 
particle during its motion increases as its initial velocity \dy(f)\ decreases. It is 
thus maximal in the limit where V^//(l) = 1 + tends towards Ve//(0m) = 
c^/fl'(0m) + g{(pm), that is in the limit g{4>m)- The time spent on the 

trajectory diverges logarithmically when approaches g{4>m) (since g' has a 
double zero at = 1). Thus, for £ 1, is exponentially close to g{(f)m)- 

We consider now the total strain constraint Eq. (lA.lip and the determination 
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of (pm- It is helpful to rewrite Eq. (lA.lip using Eq. (lA.lOp as 



S = ci + c J 



9{4>)\/Veff{(pm) - Veff{(f)) 



(A.12) 



Under this form, for large i, as 
and can be replaced by g{(pm] 
obtains 



— > g{(f)m) the integral in Eq. ( ]A.12p converges 
with an exponentially small error. Thus, one 



a/2 



+ 



9i<t>) 



(A.13) 



The existence of solutions with S £ imposes g{4>m) ^ 1 (since the integral 
term in Eq. (1A.13P is clearly positive). When g behaves as g{(j)) ~ a0'^ for 
small Eq. flA.13P reduces to 

5 ^ v^C^'^ + ^^^-'^/^ (A. 14) 



where the constant can be expressed in term of the Euler B function as 
Ccr = -8(1 — l/o", l/2)/(T. The first term in Eq. flA.14p represents a contribution 
to the total displacement that is distributed over the whole sample whereas the 
second one is a localized contribution coming from the center of the stretched 
band. If one wishes that some solutions can correspond to fractured bands, the 
second contribution should dominate the first. It should moreover be able to 
take values much larger than one (so than one can have localized solutions with 
5 3> 1). This clearly requires the exponent 1 — (j/2 to be negative and therefore 



that the function g{(f)) be chosen so that a > 2, as not ed in (IKarma et al. 



200ll ). A possible choice, made in the present work as in (IKarma et al.l . 120011 ) 



is to take a = 3 (in addition to the requirements g{l) = g'{l) = 0) and 

(,(0) = 403 - 304 (A. 15) 



For this specific choice, clearly a = 4 and C3 = 5(2/3, l/2)/3 ~ 0.862. For 
0" > 2, the r.h.s. of Eq. flA.14p has a minimum value 6c that is reached for 
4>m = 4>c with 



-i~a(fT-2)/(a(T£) (A.16) 

5, ~ y^^i^:ill£0-/2 _ (A. 17) 

a — 2 
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or more simply for our specific choice of g with a = 3, (j)c — 0.268/v^ and 
5c — 1.1 If^/"^. For an adimensionned strain 5 below 5c no non-trivial solutions 
exist. Two coincident solutions appear at 5 = 5c which separate into a stable 
lower energy solution and an unstable higher energy one when 5 > 5c a^s shown 
in Fig. lA.li The energies of the two solutions can be explicitly obtained in the 
limit i ^ 1. Eq. IA.9I can be rewritten as 

^ [Vcff{<j>J -^=J^—=2[1 - g{<f>)] (A.18) 



I \/yeff{<Pm)-Vcff{<P) 



where we have used the expression ( lA.lOp for the strip width i. For a large 
can be replaced by g{4>m) with an exponentially small error to obtain 



E 



\/2k,Sc 



d<P[l-gi<P)], 



\ 9{<f>) - 9{<Pm) 



(A.19) 



Finally, in the whole regime of interest where 5 <^ £, the phase-field minimum 
value (prn vanishes as a power of i. With the small behavior g{(j)) ~ a(j)'^, 
Eq. (1A.19|) simply reduces to 

E } 

arj - D,(t>m + 2 / d(t>[l - gm (A.20) 







where as above the constant D„ can be expressed using Euler B function 
{D„ = {a - 2)C^) and for a = 3, D3 ~ 0.862. The asymptotic form fOCTj) is 
already reasonnably accurate for £ = 3 as shown in Fig. lA.li 

As 5 becomes much larger than 5ci these two solutions correspond to the 
dominance of one of the two terms on the 1. h. s. of Eq. flA.13p . 

For the stable solution, the localized contribution to the strain dominates so 
that 

5 ^ (A.21) 



Thus, in this parameter regime, £ 3> 1 and 5 3> 5c 3> 1, 0m tends toward 
zero. As a welcome consequence, the energy of the stable solution becomes 
independent of the strain and can be identified with twice the surface energy 
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1 

= 27 = 2^2^ J d<f)^l - g{(j)) 



(A.22) 



For our specific choice of g{(f)) [Eq. (lA.lSP ] tlie numerical value of the integral 
is approximately 0.7165. As in Griffith's original theory, this stable solution 
becomes energetically favored as compared to the elastically stretched band 
when Eg becomes smaller than the purely elastic stretching energy [Eq. ( ]A.4[) ]. 



that is when A > Ap = y [27/(A + fi)]L or equivalently for 6 > 6g = 1.20^/1. 
For the unstable solution, when 6 becomes much larger than 6c, one has simply 
6 ~ V^rjH (A.23) 



The corresponding energy is simply (Eqs. flA.20p . (IA.22p ). 



2 



E^ = 2^ + J2^c -T (A.24) 



since the term proportional to 0m becomes subdominant with respect to the 
other two (and tends towards zero) when 5 moves away from 5c- In other 
terms, the energy of the unstable state is simply the energy of the elastically 
stretched band plus the energy necessary to create the the two interfaces, as 
one could have intuitively guessed. 

Finally, it is interesting to see how the profile of the fracture state depends on 
the total strain 5. The profile of a general solution is obtained from Eq. (1A.6P 

as 




(A.25) 



with as before ^ = \Jk/{28c) denotes the process zone scale. In the regime 
£ 3> 1 and 5 ^ 1 {oT equivalently 0m -C 1) the phase field profile on the 
process zone scale is independent of 5 



y 



(A.26) 



This is not true for comparable to 0m (in the regime i ^ 1 and 5^1) 
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Fig. A.l. (A) Minimum value (pm of the phase field vs. adimensionned displacement 
6 = A(/i/«;)^/^ for a strip of adimensionned half-width £ = L{2Ec/ nY^"^ = 3 (solid 
line) and in the limit of a large strip (dashed line), as given by Eq. (1A.14P with 
£ = 3. For a given displacement 5 > 6c, there are two values of (pm corresponding to 
two stationary solutions with the smallest value of (j)m for the stable one. (B) The 
energies E vs. 5 for the two branches of solutions for the same band of ^ = 3 (solid 
line), the lower branch corresponding to the stable fractured state. The asymptotic 
expression for a large band is also plotted using Eq. (|A.20p with 1 = 3 (dashed line). 
The two branches meet at a cusp, the generic behavior at a saddle-node bifurcation. 
The corresponding adimensionned energy (5^/(2^) of the (half) elastically-stretched 
band is also plotted (dotted line). It becomes larger than the energy of the stable 
fractured solution at Griffith threshold. 



where Eq. (1A.25I1 simplifies to 

y ! f p ' dp 



(A.27) 



In the same regime ~ the strain field can be written 



u 



Ll-a/2 
m 



<t>/4>n 



\ + 2p y/a 



dp 



A 



dp 



(A.28) 



Comparison of Eq. (1A.27P and Eq. (]A.28P shows that the variation of u is 
comparable to the total strain A [Eq. (]A.21] on a scale much smaller than 
the process zone length. More precisely, for different strains A, the different 
scaled strain profiles u/A are given by a unique function of y/{(t>mO- 
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B Variational equations and self-adjointness of linearized opera- 
tors 



The equilibrium phase-field equations considered in this paper are Euler- 
Lagrange equations coming from the variation of an energy density £ where 
£ depends on a set of fields ip^ (here the elastic displacements and a scalar 
phase-field) and their spatial derivatives djipai 

= o (B.1) 



We show here that quite generally for this type of equations, the allied linear 
operator is self-adjoint, a property that we used for obtaining the explicit 
formulae of section |H 

Linearization of Eq. flB.l|) around a solution ijj^^ produces the linearized op- 
erator C It is defined by its action on a set of functions Vj3 as 



d^£ d^£ [ d''£ ' 



(B.2) 



Now it is easily seen using integration by parts that for two arbitrary sets of 
differentiable functions and ly^, one has 

/d. = ji^ v^CA{M\ + boundary terms (B.3) 



The relation clearly holds separately for the first and last term on the r. h. s. 
of Eq. IB. 21 and comes from the interchange of the second and third term on 
using integration by parts. Thus, the linear operator is self-adjoint for the 
usual flat measure (here simply dx = dxidx2 on the plane). 



C Explicit computations of zero-modes and solvability integrals. 



For the convenience of the reader, we provide below some details of our com- 
putations of solvability integrals and of the Eshelby tensor line integrals. 

For plane strain, the explicit form ([T]) of the stress distribution near a crack 
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tip is conveniently obtained from the Airy function x which satisfies the bi- 
harmonic equation 



= 



(C.l) 



In polar coordinates (r, G), x is related to the strain tensor by 

CTrr = ^dl^X + ^9rX, £^60 = ^rrX ^ud CT^e = ^^OX " ^d%X (C.2) 



For a crack along the x-axis, with its tip at x = 0, the Airy function is 
determined by Eq. (IC.lll . together with zero traction boundary conditions on 
the fracture lips, arr = fre = for = ±tt. The most singular possibility 
compatible with a bounded elastic energy reads, in polar coordinates. 



r3/2 ( Xi e 36 

X = < , — [3 cos( — ) + cos( — ) 



-^[sin(-)+sin(— )] (C.3) 



The dominant divergent forms of the stress distribution follow by differentia- 
tion with the help of Eq. <^^, 
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2. + ^cos(p] (C.4) 



The relation between the strain and stress tensors (Eq. (|20|) ) and integration 
give the allied displacement field. 



Ui 



4 + K2 dl' 



r,6 



4/i V 27r 

The mode I crack tip displacement functions dl are given by 



(C.5) 



< = (5 - 8z/) cos(® ) - cos(^) 
4 = (-7 + 8z/) sin(|)+sin(^) 



(C.6) 



The corresponding mode II functions are 
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dl' = (-5 + 8u) sin(|) + 3 sin(^) 

4' = (-7 + 8u) cos(|) + 3 cos(^) (C.7) 



These expressions allow the explicit evaluation of the different integrals of 
sections H] and [5l 

i) Solvability integrals. 

For a vectorial field u = w^e^ + Meee the two components of the x-translation 
field u*^^) = dxU are 



(x) /r^N o sin(e) sin(0) 

U). ' = COS(B) OrUr OeUr H Ue 

W ^r^\ a sin(e) sin(e) 

With these formulae, one can compute the two components of the x-translation 
field u*^^''-'^ associated to the mode I the displacement field [Eq. ( ]C.5fC.6l) ]. 



^ix;i) ^ Ki _ ^^g(3Q/2) - 3 cos(e/2)] 

8/iv27rr 



^{x;i) ^ ^ sin(3e/2) + 3 sin(e/2)] (C.9) 

8/iv27rr 

The corresponding strain tensor reads 



= drv^r'^ = - ' [(7 - 8Z.) cos(3e/2) - 3 cos(e/2)] 

lo/iv27rr'^'^ 

4t^= ^-[u^:-^'^ + deut''^] =----|L-^[(1-8z/)cos(30/2) + 3cos(0/2)] 



u 



ix;I) 
rB 





16/i 














3Ki 


16/i 





lidJe''^ + deu^;-''^] = -^^^[sin(3e/2) + sin(e/2)] (C.IO) 



2' 

The two needed components of the allied stress tensor follows 
a^'^ = p [-7cos(3e/2) + 3cos(G/2)] 

aj^/) = ^\,, [sin(3e/2) + sin(e/2)] (C.ll) 
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With these exphcit expressions, one can evaluate the solvabihty integrals for 
the perturbed displacement field ui\ 



u 



(1) 



4/i V 27r 



5K^ 4 + 6K2 4' 



i = r,Q 



(C.12) 



With Eq. flC.4yC.9p and flC.6yC.lll) . one obtains for the two integrals, 



Ki5K 



8/i 



i[5-6z/] 



8// 



where a^j^j^ is the perturbation of the stress tensor corresponding to ul''' . It is 
given by the mode I part of Eq. flC.4p with Ki and K2 replaced by 5Ki and 



(C.13) 



(1) 



5K2. For symmetry reasons, only the mode I part of and a^j' contribute 



1-3 



to the integrals. Substracting the two Eqs. fIC.lSp . one finally obtains 



■[l-u] 



(C.14) 



which is equation fl27P of the main text. 

The corresponding expressions for a y-translation field u*^^-* = dyU associated 
to a vectorial field u = w^e^ + uqCq are 



u 



(u) ■ n C0S(9) ^ COS(O) 

^y> -= sm(e) drUr H — deUr — ue 



(v) ■ n cos(6) „ cos(6) 
Uq = sm(B) OrUQ H OeUe H Ur 



(C.15) 



This gives for the two components of the y-translation field u*^^'-'^^ associated 
to the mode I displacement field of Eq. flC.6P 



u 



u 



e 



Tcr 



with the corresponding strain tensor 



[(7 - 8z/) sin(3e/2) - sin(e/2)] 
[(5 - 8u) cos(30/2) - cos(0/2)] 



(C.16) 



u 



16/iv^r3/2 



[(7 - 8z/) sin(3e/2) - sin(e/2)] 
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= - 16^^,3/2 - sm(3e/2) + sin(e/2)] 



4^/) = 77— ^17^[3 cos(3e/2) + cos(e/2)] (C.17) 
and stress tensor 



a^'^ = [-7sin(3e/2) + sin(e/2)] 

a?/) = [3cos(3e/2) + cos(e/2)] (C.18) 



This gives for the two integrals of interest 



<i»."l'V|f" = :?M2|-3 + 2.] (C,19) 



where only a^f and u'l^\ the mode 11 parts of the perturbed the stress tensor 
alf and displacements fields u\^^ (Eq. (1C.12P contribute to the integrals. 



One finally obtains for the difference of the two integrals of Eq. (]C.19|) 



ci.,[.rV«-n«o-r^] = ^[l-.] (C.20) 



which is equation fl30|) of the main text. 

ii) Components of the configurational force on the crack tip 

The two line integrals [Eq. ( I38|) ] giving the two components of the configu- 
rational force on the crack tip can be directly evaluated with the help of the 
above results. 

y^con/) ^ j^^jE^^^ j ds [Estrainni- CTjkuf] (C.21) 
A->B A~*B 



where u*^') = diU is the translation field in the direction i{i = x,y), with mode 
I and mode II components included 

= 4*'^) + 4*'^^) (C.22) 
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The mode I components m^*'^^ are given by Eq. (\C.9\i and flC.16p . A similar 
computation gives their mode II components with the help of Eqs. flC.5IIC.7l) 
and <^M> 



4-'''^ = — ^ [(-7 + 8u) sin(3e/2) + sin(e/2)] 
8/iv27rr 

{x;II) _ K2 



^y.,,,, ^ z ^ cos(3e/2) + cos(e/2)] (C.23) 

8/iv27rr 



and 



^{yji) = [(7 - 8z/) cos(3e/2) + 5 cos(e/2)] 

8/iv27rr 

^{y,ii) ^ K2 ^ sin(3e/2) - 5 sin(e/2)] (C.24) 

8/iv27rr 



With these formulae and the stress tensor expression Eq. flC.4p . one obtains 



+7r 

Jrde [arrui''^ + a,eM^"^] = —[Kf{-3 + 2v) + Kl{-b + 6z/)] 

— TT 

yrc/e [a^M^^) + areui>] = — 7^1^2(3 - 2u) (C.25) 



Furthermore, with Estr-ain = {crijaij — uauajj) and the stress tensor expression 
Eq. (]C.4p . one obtains 



+7r 

/rrfe cos{Q)£,tra^n = —{Kl-Kl){l-2v) (C.26) 
J 8/i 

-TT 

yrde sin(e)£,i,,„ = -— i^iir2(l - 2z/) (C.27) 



Substraction of Eq. (1C.25I) from Eq. ( 1C.26I) gives the usual expression of 

p(con/) 

picnf) ^ [cos(e)^,,,,,„ - a„u^:^ - a^eu^^^] = ^(^1 + i^l)(C.28) 

— n 

as given in Eq. (HP]) in the main text. 
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The other component of the configurational force F2^°"'^^ is similarly obtained 
by substracting Eq. flC.25D from Eq. flC.27p . 



1 

piconf)^ jdQ[sm{Q)£^^^,^^^-arrU^y^ -areut\ = ^^1^2 (C.29) 



which is Eq. (jUj) of the main text. 
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